function sandbox = fsb_HRF_auto_adjust(idat,sandbox,param)

% FSB - DEV: shift predictors in such a way that the correlation with
% the chosen ROI is maximal
%
% EXAMPLE:
% sandbox = fsb_HRF_auto_adjust(idat,sandbox,param)
%
% INPUT:
% idat:             4D image data
% sandbox:          sandbox experiment struct
% param:            parameter containing information on the ROI and chosen predictor
% param.slice_n:    Position of the ROI
% param.roi:        Size of the ROI
% param.hrf_pred:   Hemodynamic predictor chosen
%
% OUTPUT:
% sandbox:          adjusted sandbox experiment struct
% 
% CALLED BY:
% FSB.m
%
% NOTES:
% in development, use at your own risk.
% This function is not to be used to improve statistical results in your
% area of interest
%
% Copyright 2010 MPI for Biological Cybernetics
% Author: Steffen Stoewer
% License:GNU GPL, no express or implied warranties
% 
% $ Revision 0.1
%
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

if isfield('fullscan','sandbox') && sandbox.fullscan == 0
    disp('Concatenated Trial, aborting...')
    return
end

h = waitbar(0,'Finding optimal shift...');
testfac = zeros(100,2);
onset_spare = sandbox.onset;
intrial_spare = sandbox.intrial;
hemodynamics_spare = sandbox.hemodynamics;
inparg = struct;
dim = size(idat);

%~~~~~~~~~~~~~~~~~~~~~~~
% Check for old version incompatibilities
%~~~~~~~~~~~~~~~~~~~~~~~

try
    sandbox.volmat.evt_onsets = sandbox.origdata.volmat.evt_onsets;
    sandbox.intrial = sandbox.origdata.intrial;
    sandbox.hemodynamics = sandbox.origdata.hemodynamics;
    sandbox.timemat.evt_onsets = sandbox.origdata.timemat.evt_onsets;
    sandbox.onset = sandbox.origdata.onset;
catch
    disp('No prior optimization found')
end

try
    inparg.evt_dur = sandbox.volmat.evt_dur;
    inparg.evt_onsets = sandbox.volmat.evt_onsets;
catch
    try
    sandbox = fsb_volmat(sandbox);
    inparg = sandbox.volmat;
    catch
        disp('sandbox.volmat not defined yet, please set up experiment struct first')
        return
    end
end

inparg.TR = sandbox.TR;
inparg.volumes = dim(4);
inparg.deriv = 0;

%~~~~~~~~~~~~~~~~~~~~~~~
% Shift experiment information back and forth and look for maximal signal
%~~~~~~~~~~~~~~~~~~~~~~~
for x = 1:100
    sandbox.intrial = intrial_spare;
    sandbox.onset = onset_spare;
    sandbox.hemodynamics = hemodynamics_spare;
    y = x/10-5;
    waitbar(x/100);
    shiftparam = y;
    if shiftparam>cellfun(@min,sandbox.volmat.evt_onsets)
        for y=1:size(sandbox.volmat.evt_onsets,2);
            inparg.evt_onsets{y,1} = sandbox.volmat.evt_onsets{y}+shiftparam;
        end

        [U,bf,p,sandbox.hemodynamics] = fsb_conv_onsets_time(inparg);
        sandbox.intrial = circshift(sandbox.intrial,[round(shiftparam) 0]);

        trial_idat = idat(:,:,:,sandbox.intrial(:,1)>0);
        trial_pred = sandbox.hemodynamics(sandbox.intrial(:,1)>0,param.hrf_pred);

        try
            trial_hrf = squeeze(mean(mean(mean(trial_idat(param.slice_n(1)-param.roi:param.slice_n(1)+param.roi,...
                param.slice_n(2)-param.roi:param.slice_n(2)+param.roi,param.slice_n(3)-param.roi:param.slice_n(3)+param.roi,:)))));
        catch
            errmsg = lasterror;
            disp('Error in file:'); disp(errmsg.stack(1).file); disp('Line:');disp(errmsg.stack(1).line);
            disp('ROI edge exceeds matrix dimensions, aborting operation')
            return;
        end

        testfac(x,1) = corr(trial_hrf,trial_pred);
        testfac(x,2) = y;
    end

end

%~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Now restore everything to the original state
%~~~~~~~~~~~~~~~~~~~~~~~~~~~

sandbox.intrial = intrial_spare;
sandbox.onset = onset_spare;
sandbox.hemodynamics = hemodynamics_spare;

testmax = max(testfac(:,1));
testfin = find(testfac(:,1)==testmax);
try
    shiftparam = testfin/10-5;
    disp(['Optimal shiftparameter ' num2str(shiftparam) ' Volumes' ]);
catch
    disp('Shift not useful, no improvement in data');
    close(h);
    return;
end;

close(h);
%~~~~~~~~~~~~~~~~~~~~~~
% Apply optimal parameter to data
%~~~~~~~~~~~~~~~~~~~~~~
sandbox.origdata.volmat.evt_onsets = sandbox.volmat.evt_onsets;
sandbox.origdata.intrial = sandbox.intrial;
sandbox.origdata.hemodynamics = sandbox.hemodynamics;
sandbox.origdata.timemat.evt_onsets = sandbox.timemat.evt_onsets;
sandbox.origdata.onset = sandbox.onset;
for x=1:size(sandbox.volmat.evt_onsets,2);
    inparg.evt_onsets{x,1} = sandbox.volmat.evt_onsets{x}+shiftparam;
    sandbox.volmat.evt_onsets{x,1} = sandbox.volmat.evt_onsets{x}+shiftparam;
end

[U,bf,p,sandbox.hemodynamics] = fsb_conv_onsets_time(inparg);
sandbox.intrial = circshift(sandbox.intrial,[round(shiftparam) 0]);
disp('Intrial shifting done');
sandbox.onset = circshift(sandbox.onset,[round(shiftparam) 0]);
disp('Onset shifting done');

end
